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Abstract 

In the presented paper a family of angiogenesis models, that is a generalisation of the Hahn¬ 
feldt et al. model is proposed. Considered family of models consists of two differential equations 
with distributed time delays. The global existence and the uniqueness of the solutions are proved. 
Moreover, the stability of the unique positive steady state is examined in the case of Erlang and 
piecewise linear delay distributions. Theorems guaranteeing the existence of stability switches 
and occurrence of Hopf bifurcations are proved. Theoretical results are illustrated by numerical 
analysis performed for parameters estimated by Hahnfeldt et al. {Cancer Res., 1999). 
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1 Introduction 


Angiogenesis is a very complex process which accompanies us through our whole life starting from 
the development of the embryo and ending with wound healing in an advanced age since it is a pro¬ 
cess of blood vessels formation from pre-existing structures. Clearly, it is often considered as a vital 
process involved in organism growth and development. On the other hand, it can also be a patho¬ 
logical process since it may promote the growth of solid tumours. Indeed, in the first stage of solid 
tumour development cells create multicellular spheroid (MCS) — a small spherical aggregation. For 
a fast growing tumour (comparing to the healthy tissue) tumour’s cells located in the centre of MCS 
receive less and less nutrients (such as glucose and oxygen) since they are only supplied through 
the diffusion of that substances from the external vessels. Hence, when MSC reaches a certain size 
(usually around 2-3 mm in diameter, |T]|2j) two processes are observed: the saturation of growing cel¬ 
lular mass, and the formation of necrotic core in the centre of MCS. The tumour cells that are poorly 
nourished secrete number of angiogenic factors, e.g. FGF, VEGF, VEGFR, Angl and Ang2, which 
promote the proliferation and the differentiation of endothelial cells, smooth muscle cells, fibroblasts 
and also stabilise new created vessels initiating the process of angiogenesis. From that point of view, 
angiogenesis can be also considered as an essential step in the tumours transition from less harm for 
hosts avascular forms to cancers that are able to metastasise and finally cause the lethal outcome of 
the disease. 

In [TJ Folkman showed that the growth of solid tumours strongly depends on the amount of blood 
vessels that are induced to grow by tumours. He summarised that, if the tumour could be stopped 
from its own blood supply, it would wither and die, and nowadays he is considered as a father of an 
anti-angiogenic therapy the approach that aims to prevent the tumour from developing its own blood 
supply system. 

On the other hand, angiogenesis has also an essential role in the tumours treatment during chemother¬ 
apy since anti-cancer drugs distributed with blood need to be delivered to the tumour and efficiently 
working vessels net allows the drugs to penetrate better the tumour structure. One also needs to take 
into account the fact that the vessel’s net developed due to angiogenesis initiated by tumour cells is 
not as efficient as the net in healthy tissue e.g. consists of loops. That causes large difficulties during 
the treatment of tumours because drugs transport is often ineffective. Because of the reasons ex¬ 
plained above, the process of angiogenesis is very important for the solid tumour growth and also for 
the anti-tumour treatment. Hence, a number of authors described that process using various math¬ 
ematical models: macroscopic mu- individual-based models [ fl2[|T3| or hybrid models |14j|T5j. 
One of the most well recognised angiogenesis models was proposed by Hahnfeldt et al. [JTTJ later on 
studied in [8]. Among others, models with discrete delays describing the same angiogenesis process 
were recently considered in [5j and [16]. Models investigating the tumour development and angio¬ 
genesis in the context of anti-angiogenic therapy, chemotherapy or radiotherapy were also considered 


in 117 -231, in the context of optimal treatment schedules [24]-29] or in context of vessels maturation 


in [30 311. 

Let consider the family of the angiogenesis models that consists of two differential equations in 
the following form 


l p(,) = rp(t)h [*^ 


(1.1a) 

(1.1b) 


where variables p(t) and q(t) represent the tumour volume at time t, and the maximal tumour size for 
which tumour cells can be nourished by the present vasculature. Thus, the ratio ^ is interpreted as 
the measure of vessels density. The function h reflects dependence of tumour growth rate on the vas- 
cularisation. We assume h is decreasing. In the literature, in the context of model (| 1.1 1), it is usually 
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assumed that the tumour growth is governed by logistic, generalised logistic (see |[8|-[10|) or by Gom- 
pertzian (see [8}|TTJ) law. The parameter r reflects the maximal tumour growth rate. The delays T\ 
and T 2 present in the system represent time lags in the processes of the tumour growth and the vessels 
formation, respectively. Moreover, it is assumed that the population of the endothelial cells depends 
on both: the stimulation process initiated by poorly nourished tumour cells and the inhibiting factors 
secreted by tumour cells causing the lost of vessels. The parameter a reflects the strength of the de¬ 
pendence of the vessel dynamics on the ratio Parameters ciu and b are proportionality parameters 
of inhibition and stimulation of the angiogenesis process, respectively. We consider two processes, of 
which one appears on the surface of a sphere and the other the inside the sphere, and due to the spher¬ 
ical symmetry assumption we suppose that the ratio of the surface process depends proportionally to 

2 

the tumour surface and it is represented by the exponent 2/3. Thus, the term -a H q(t)pi(t) appears in 
Eq. ( |l.lb[ ). In model ( |1.1[ ) the spontaneous loss of functional vasculature is represented by the term 
-pq(t). In [IT], the parameter p was estimated to be equal to zero. On the other hand, the term -pq(t) 
for p > 0 can be also interpreted as a constant continuous anti-angiogenic treatment, hence it is also 
considered in the presented study. For detailed derivation of system ( |1.1[ ) without delays and with 
logarithmic function h from a reaction-diffusion system see CD- 

Model ( | 1.1 1) with Gompertzian tumour growth, a = 1, and without time delays was firstly pro¬ 
posed by Hahnfeldt et al. in HD- The modification of the Hahnfeldt et al. model was considered 
by various authors. d’Onofrio and Gandolfi [8] proposed system of ordinary differential equations 
(ODEs) based on Hahnfeldt et al. idea but with h being linear or logarithmic function, and with 
a = 0, that is when the dynamics of the second variable of the model does not depend on the vas- 
cularisation of the tumour. Later in [32] Bodnar and Forys introduced discrete time delays into the 
model with the Gompertzian tumour growth and the family of models for parameter a from the inter¬ 
val [0,1] was considered. The analysis was later extended in [171. In the same time, d’Onofrio and 
Gandolfi in [10] studied the model for a = 0 with discrete time delays and different functions h. The 
analysis of the family of the models with discrete delays and the Gompertzian tumour growth was 
extended in [33] 341 where the directions of the existing Hopf bifurcations were studied. 

However, in all these papers only discrete time delays were considered. One of the reasons is the 
fact that mathematical analysis of such models is easier than the models with distributed delays. Of 
course, in reality the duration of a process is never constant and it usually fluctuates around some 
value. Thus, we believe that delay is somehow distributed around some average value. This is the 
main reason why in this paper we study the modification of the family of models ( |1.1[ ) in which the 
delays are distributed around theirs average values r { and r 2 instead of being concentrated in these 
points. On the other hand, we should also point out, that it might be difficult to estimate a particular 
delay distribution form experimental data. 

Up to our knowledge such modification of the family of models has not been considered yet. We 
are also not aware of any result considering linear systems similar to the linearised system of the 
model with distributed delays considered in this paper. In fact known analytical results regarding 
the stability of trivial steady state for linear equations with distributed delays are rather limited, see 
e.g. [135] and references therein. Most results concern a single equation (see [[36|-38|, and references 
therein). In some papers a second order equation or a system of equations with distributed delays are 
considered (see [ 391 for the second order linear equation, and [001 for the Lotka-Volterra system). 
However, in these two cases the time delays are only finite and the delayed terms do not appear on 
the diagonal of the stability matrix, and hence these results cannot be applied directly to the system 
considered in this paper. The single infinite distributed time delay is considered in [41 ] for the model 
of immune system-tumour interactions, however in that paper only an exponential distribution is 
considered and the linear chain trick that reduces the system with infinite delay to a larger system of 
ODEs is used. 

A part of our results is based on application of the Mikhailov Criterion, which generalised version 
is formulated and proved in the Appendix. This generalisation plays an essential role in studying local 
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stability of the steady state with certain distributed delays. 

The paper is organised as follows: in Section[2]considered family of angiogenesis models with dis¬ 
tributed delays is proposed, and a proper phase space and initial conditions are defined. In Section [3] 
the mathematical analysis including the global existence and uniqueness of solutions, the stability of 
the existing steady state of considered family of models for different types of considered delay distri¬ 
butions is presented. In that section we also discuss possibilities of stability changes of the positive 
steady state. Next, our stability results are illustrated and extended by numerical simulations. In Sec¬ 
tion [5] we discuss and summary our results. Finally, for completeness in [A] we formulate generalized 
Mikhailov Criterion used in the paper and we prove it. 


2 Family of angiogenesis models with distributed delays 


The effect of the tumour stimulation of the vessels growth as well as the positive influence of the 
vessel network on the tumour dynamics is neither instantaneous nor delayed by some constant value. 
Hence, to reflect reality better, instead of constant delays considered earlier in j8], [ 331 or [34] we 
consider the delays that are distributed around some mean value and their distributions are given by 
the general probability densities /). We study the following system 


pit) = rp(t)h (fi (r) —— dr |, 


dr 

_d 

dr 


q{t - t ) 

q(t) = <?(r)| b f 2 (r) (^^ ) dr - a H p 2/3 (t) - p |, 


(2.1a) 

(2.1b) 


,q(t - t ) 

where fj(s) : [0, oo) — > R> are delay distributions with the following properties: 

J f~*oo 

Ms)ds = 1, i = 1,2 and 

o 

sfi(s) ds < oo, i = 1,2. 
o 

Moreover, we assume that the function h has the following properties: 

(PI) h : (0, oo) —» R is a continuously differentiable decreasing function; 

(P2) h( 1) = 0; 

(P3) h\ 1) = -1. 

Note, that properties (P2) and (P3) do not make our studies less general as a proper rescaling and 
a suitable choice of parameter r always allows us to arrive to the case /?(1) = 0 and h'( 1) = -1. 
Note also, that an arbitrary probability distribution defined on [0, +oo) with a finite expectation fulfils 
assumptions (HI) and (H2). 


To close model ( |2. 1 1 ) we need to define initial conditions. We consider a continuous initial func¬ 
tion (p : (— oo, 0] —> R 2 . For infinite delays we need to regulate the behaviour of this function as t tends 
to -oo. To this end, we introduce a suitable space. Let us denote as C - C((-oo, 0], R 2 ) the space 
of continuous functions defined on the interval (-oo,0] with values in R 2 . For an arbitrary chosen 


non-decreasing continuous function q : (-oo, 0] 
space 


R + such that lim q(0) = 0, we define the Banach 


K n 


(fi e C : lim (f>i6)qiff) 

0 —^—oo 


0 and sup \(fii6)qid)\ < oo 

0e(—oo,0] 
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with a norm 

IMIi/ = sup \(fi(9)r](d)\, for any <p e <K n . 

0e(—oo,0] 

We need initial functions 0 to be in f K n for some arbitrary non-decreasing continuous function r\ that 
tends to 0 for arguments tending to -oo (see [[42 ]). If the delay is finite, that is if supports of f\ and f 2 
are compact, we suppose initial functions 0 e C([-r max , 0], R 2 ), where the interval [-r max , 0] contains 
supports of f, / = 1,2. However, this is equivalent to consideration of the space 'K n with an arbitrary 
tj being equal to 1 on [-r max , 0] and decreasing to 0 in -oo. If the support of the probability density 
is unbounded and the initial function (p is also unbounded, than we need to choose an appropriate 
function r\ (which in turn implies a choice of the phase space ’K). The function // must be chosen 
to control the behaviour of the initial function in -oo. However, due to biological interpretation of 
the model parameter, it is reasonable to restrict our analysis to the globally bounded initial functions. 
Thus, we can chose an arbitrary positive continuous non-decreasing function q such that t](0) —> 0 
for 6 —> -oo for example // = e°. Nevertheless, because such assumption would not make theorems’ 
formulation simpler, we decided to present and prove theorems for arbitrary initial functions and not 
to assume a particular form of the function q. 


3 Mathematical analysis of family of angiogenesis models with 
distributed delays 


Due to properties (P2) and (HI) the steady states of system (2.1) are the same as for Hahnfeldt el 


ol. and d’Onofrio-Gandolfi models without delays or with discrete delays. Thus, system ( |2.1[ ) has 

a unique positive steady state (p e ,q e ), where p e - q e - (^) 2 (compare e.g. [331) if and only if b > p. 
Hence, in the rest of the paper we assume that b > p holds. 


Note, that solution to the equations defined by ( |2.1[ ) can be written in the exponential form. Thus, 
using the same argument as in [ 101 or in [33 [ (for the discrete delay cases), we deduce that R 2 is the 


invariant set for system ( |2.1[ ). 

Note, that in system ( |2.1| ) terms with delay are of the form p/q. This, together with the invariance 
of R 2 suggests the following change of variables 


In | — 


In 


Me 

qPe 


giving 


^X(t) = rh /i(t) e- v( ' r) drj , 

( s~*CO \ /~\00 

J f] (t) e y(, ~ T> drj - b J Mr) d ry< ‘~ T> dr + (b - p) e^ x(,} +p. 


d 

dr 


(3.1) 


Newly introduced variable y allows us to consider the system where only one variable (y) has delayed 
argument, whereas in system (|2.1[) both variables have delayed arguments. Clearly, the steady state 


for re-scaled system (3.1) is (x e ,y e ) = (0,0). 


It should be mentioned here that the scaling procedure transforms R- into the whole R~, so space 


'K n is the appropriate phase space for the model ( |3. 1 [ ). 


3.1 Existence and uniqueness 

Proposition 3.1 Letq : (—oo,0] —» R + be a continuous, non-decreasing function such that lim t](0) = 

9 —>—oo 

0, let functions f fulfil (H1)-(H2), and let function h fulfils (P1)-(P3). For an arbitrary initial func¬ 
tion (p = (0i,02) 6 'K n there exists t# > 0 such that system ED with initial condition x(t) = 0|(/), 
y(t) = 0 2 (t)/or t 6 (-oo, 0], has a unique solution in 7C,, defined on t e [0, tf). 
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Proof : The right-hand side of system ( |3.1[ ) fulfils a local Lipschitz condition. In fact, the assumption 
(PI) implies that the derivative of h is bounded on an arbitrarily chosen compact set of (0, +oo). Thus, 
the function h is locally Lipchitz continuous function as well as all functions on the right-hand side 
of p.l| ). This implies that the Lipschitz condition is fulfilled for any bounded set in 'K n . Hence, there 
exits a unique solution to system ( |3.1[ ) defined on t £ [0, tf) (see [|42j Chapter 2, Theorem 1.2]). ■ 


Theorem 3.2 (global existence) If assumptions of Proposition [XT] are fulfilled, the probability den¬ 
sities f and an initial function are globally bounded, then solutions to 63' are defined for all t > 0 . 

Proof : The right-hand side of ( |3 .1 1 ) can be written in a functional form as 

d d 

—x(t) = G\(x t ,y t ), —y(t) = G 2 (x t ,y t ), 

at at 

with 

Gi (<f) x , (p y ) = rh fi (t) c (M ~ TI drj , 

( z>oo \ Z"»00 

J o f\ (t) e ,/ ’ v(_T) drj - d J o Mr)a a ^-r) dT + ( b-p) e l^+p, 

for (fix, <p y ) 6 C. Note that due to (PI) we have \G\(cp x , </ v )l < rh(exp(-\\f y \\)) and similar inequality 
holds for | G 2 (f x , f y )\- This implies, that the function (G\, G 2 ) maps bounded sets of C into bounded 
sets of R 2 (so their closures are compact). Thus, if the solution to (|3. 1 [) cannot be prolonged beyond 


the interval [0, T ), for some finite T, then lim(||x,|L + ||y,||, ; ) = +oo (see [42 Chapter 2, Theorem 2.7]). 

t—>T 

In the following, we show that x and y are bounded on [0, T) hence, the solution exists for all / > 0. 

J ~oo 

J)(t) dr > 0 for i = 1,2. Using the step method we 
choose the time step equal to 8 and we show that the solution to ( |3.1[ ) can be prolonged on the interval 
[0, <5] and hence on the interval [n8, (n + 1 )8 ]. for any n 6 N. Let letters C, denote constants that will 
be chosen later in a suitable way. 

First, we provide an upper bound of the solution. Due to boundedness of y{t) for t < 0, we write 

/"*oo r* oo r* oo 

I Mr) e - v "" T1 dr = I /,(r) e y(f_r) dr + I ./Hr) e v, '" r ' dr > I /Hr) e- v "‘ n dr > C,, 

Jo Jo Js Js 

for all t £ [0, d]. This estimation leads to a conclusion 

d 


dr 


x(t) < rh(C i) 


x{t)<C 2 forre[0, d]. 


From the second equation of (3.1), arguing in a similar way, we have 


—y(r) < rh{C\) + (b - p) M +p 


y(t) < C 4 , forte [0,d]. 


Now, we proceed with a lower bound of the solutions. Splitting integral on the interval (0, +oo) into 
two integrals and using the fact that y(t) is bounded we obtain 

J r*oo /~*S s~*co s\S 

f(r) Q y{t ~ T) dr = /Hr) e* f - T) dr + /Hr) dr < /,(r) e Vw dr + C, 

o Jo Js Jo 

t 

f (r) e >M dr + C 5 < e VM +C 5 = C 6 , 


J r»o 

0 


where yM is an upper bound of y{t) on the interval [-r, d - r]. This estimation together with a similar 
argument applied to the second integral in the second equation of (|3.1[) yields 


0 


Mr) e 


aytt-T) 


dr < C 7 . 
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Therefore, from the second equation of (|3.1[) we have 


-y(t)>rh(C 6 )-bC 7 => y(f)>C 8 , for te[0,6]. 

The boundedness of y together with the form of first equation of ( |3.1[ ) implies boundedness of x(t) on 
[0, d]. Now, knowing that x(t) and y(t) are bounded on the interval (- 00 , d], an analogous reasoning 
allows us to deduce the boundedness of the solution on the interval (- 00 ,2d]. Hence, the mathemat¬ 


ical induction yields the boundedness of the solutions to ( |3.1[ ) on each compact interval included in 
[ 0 , + 00 ), which completes the proof. ■ 


3.2 Stability and Hopf bifurcations 

In this section we study the local stability of the steady state (0,0) using standard linearisation tech¬ 
nique. Linearisation of system (|3.1[) around the steady state (0,0) has the following form 


df 


x(t) 


/~*o 

= rh\ 1 ) 

Jo 


f\ (T)y(t - t) dr. 


d 2 r°° 

y(t ) = ^(b - p)x(t) + J (r/T(l)/i(r) - abf 2 (r))y(t - r) d t, 


(3.2) 


and, due to the equality h'( 1) = -1 (see Property (P3)), the corresponding characteristic function is 
given by 

/~*oo 

J r ( /i(r) e~' lr dr 

Jo 

r \ /~*00 

-^(b -p) A + J (r/i(r) + abf 2 (r))e~ AT dr 


W(d) = det 


Thus, 


W(d) 


A 2 + A 


(rf\ (r) + abf 2 {T))e~ AT dr + y (b - p) J /, (r) e _/lT dr. 


In general, as the probability densities one considers the specific distributions which describe the 
experimental data or studied phenomena in the best way. In the present paper, we consider two 
particular types of the probability densities. One type is given by 


1 


cr + r, r e [cr - s, cr), 


f(r) = — <s + cr - r, T 6 [cr, cr + e], 


for cr > s,i = 1 , 2 , 


(3.3) 


0 otherwise. 

We call probability densities defined by ( |3.3[ ) piecewise linear probability densities. Note, that for 


s —» 0 the functions f converge to Dirac delta at the points cr and therefore, system ( |2.1[ ) becomes 
a system with discrete delays considered by [j8j[T7 33 341. Since we assume that all considered 
probability densities are defined on interval [0, 00 ) condition cr > s must be fulfilled. Note, that for 
cr > £ the average value of f is equal to cr and the standard deviation is equal to e/ V6. For cr < £ 
one obtains so-called neutral equations, for details see [43). 

In this paper, we also consider second type of probability densities so-called the Erlang probability 


densities separated or not from zero by cr > 0 , i.e. we study system ( | 2 . 1 | ) with functions given by 

Mr) = g m {r ~ cr ), / 2 (r) = g m {r - cr ), (3.4) 

where g m ,(r), i = 1,2, are called non-shifted Erlang distributions and are defined by 


g mi (s) = 


( m i - 1 )! 


(as) 


nil—l -as 

e , 


(3.5) 
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with a > 0, s > 0. The mean value of the non-shifted Erlang distribution g nh is given by —, while 
the variance is equal to Hence, the average delay is equal to this mean and the deviation 
measures the degree of concentration of the delays about the average delay. Clearly, the non-shifted 
Erlang distribution is a special case of the Gamma distribution, where the shape parameter m, is 
an integer. It is also easy to see that the non-shifted Erlang distribution is the generalisation of the 
exponential distribution, which one gets taking m, = 1. On the other hand, for the non-shifted Erlang 
distributions when m, —» +oo the probability densities g m converge to a Dirac distributions and hence 
the system with discrete delays ( |1.1| ) is recovered as a limit of Erlang distributions for system ( |2. 1 [ ). 
For the shifted Erlang distribution the mean value is cr + -, while variances stay the same as for the 
non-shifted case. 


3.2.1 Erlang probability densities 

The Erlang probability densities separated from zero by cr read 


fi(r) 


a m '(T - &) m ‘- x 


-a(r-o-) 


, for t > cr 


0 n i - 1)! 

and 0 otherwise, for i = 1,2. Note, that in this paper we consider only the case a > 0. Then for 
probability densities given by ( |3.4| )-( |33] ) we have 

f(j) e“' ir dr 


f 


a ‘ 


-Acr 


(a + A) n 

Hence, the characteristic function has the following form 


TO) = A Z +Air 


mi 


-Act 


7 m 2 


(a + A)' n ' 


+ab- 


-Act 


0 a + A)’" 2 


2 r 

+ -(b-p) 


7 mi 


-Act 


(a + A) m i 


Define 


(3 = r + ab and y 


2 r(b - /i) 


(3.6) 


(3.7) 


Proposition 3.3 Let b > /a. The trivial steady state of system (j3TJ) with the non-shifted Erlang 
probability density given by ( |3.4[ )-( [33| ) (cr = 0) is locally asymptotically stable if 


(i) a(3 > y for m { = m 2 = 1 ; 

(ii) a > for ni\ = m 2 = 2 ; 

(iii) a > |/3 and 8/3a 3 - 3(8y + 3f3 2 )a 2 + 3yf3a - y 2 > 0 for m\ = m 2 = 3; 

(iv) 2 a(a + r) > (a(ci + ab) + y)+ (a( a+X+ y) f or = 1 and m 2 = 2 ; 

(v) a > f/3 + jr - ab for m\ = 2 and m 2 = 1 . 

Proof : For cr = 0 we clearly have 


TO) = A 2 + Air 


a m , a”' 2 \ 2 r , a m ' 

-1- ab - H- (b - u.) - 

(a + A) mi (a + A) m2 ) 3 (a + A)" !l 


As it can be seen, the advantage of using the Erlang distributions (not separated from zero) is that 
instead of studying existence of the zeros of the characteristic function one can study existence of 
roots of a polynomial, thus the stability analysis is easier. 
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First, consider m\ - m 2 = m. We investigate the behaviour of the roots of polynomial 

A 2 (a + A) m + a m (A/3 + y). 


(3.8) 


where/3 and y are given by ( |3.7[ ). 

The Routh-Hurwitz stability criterion (see eg. p4| ) gives that the necessary and sufficient condi¬ 
tions for the stability of trivial steady state of system ( |3.2[ ). Clearly, the degree of polynomial ( |3.8[ ) is 
larger for larger m. However, in each case of equals m’s (sometimes tedious) algebraic calculations 
lead to some conditions that guarantee the stability of the considered steady state. 

For m = 1 (i.e. for the exponential distribution not separated from zero) we have a polynomial of 
degree three 


A~(a + A) + a (A/3 + y), 


while for m = 2 we have 


T 4 + 2 aA i + a 2 A 2 + fia 1 A + ycr . 


The case m = 3, is the most complicated since a direct calculation shows that ( ]3.8| ) is the polynomial 
of degree 5 

A 5 + 3a A 4 + 3 a 2 A 3 + a 3 A 2 + a s (3A + a 3 y. 

For m i = 1 and m 2 = 2, polynomial ( |3.8[ ) reads 

A 4 + 2 a A 3 + a(a + r) A 2 + [a 2 f5 + yaj A + ya 2 , 

while for mi = 2 and m 2 = 1 we have 

A + 'l.ci A + ala + ab) A + cr[5 A + ycr . 


Proposition 3.4 For b > p, cr = 0, m { = 1, m 2 = 2 and 


ab 2 y 

r < ab or r > — H—-, 
2 ab 


(3.9) 


there exists a > 0 such that for a > a the trivial steady state of system ( |3.1[ ) with the non-shifted 


Erlang probability densities given by ( |3.4[ )-( |3.5[ ) is locally asymptotically stable and it is unstable for 
a G (0, a). 


Proof : Since a{a + ab) + y > 0 holds, the stability condition for the case mi = 1, m 2 = 2 (the 
condition (ii) of Theorem [33] ) is equivalent to 

W a (a ) = a 4 + 2r a 3 + [abllr - ab) - 4 y) a 2 + 2 y(r - ab) a - y 2 > 0. 


The Descart’s rule of signs implies that if one of conditions (3.9) holds, then there exists exactly one 


simple positive real root of W a (a). We denote it by a. Then the trivial steady state system ( |3.1[ ) is 
locally asymptotically for a > a and unstable for 0 < a < ci. ■ 


Theorem 3.5 Considered system (3.1) with the shifted Erlang probability densities given by (3.4)— 


(3.5). 


(a) Let nij = m, i = 1,2. If the trivial steady state is 

(i) unstable for cr = 0 , then it is unstable for any cr > 0 , 
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(ii) locally asymptotically stable for cr = 0, then there exits cr 0 > 0 such that it is locally stable 
for all cr 6 [0, cr 0 ) and unstable for cr > cr 0 . At cr - cr 0 the Hopf bifurcation occurs. 

(b) Let mi = 1, m 2 = 2. If the trivial steady state is locally asymptotically stable for cr - 0, then 
there exists exactly one cr 0 > 0 such that for cr e [0, <x 0 ) it is locally asymptotically stable and it 
is unstable for cr > cr 0 . At cr = cr 0 the Hopf bifurcation occurs. 


(c) Let nii = 2, m 2 = 1. 

(i) If the steady state is locally asymptotically stable for cr = 0, and the function 

F{u) = u 4 + 2 a 2 u 3 + u 2 a 2 (a 2 - a 2 b 2 ) - ua 2 [af> 2 - 2abyj - y 2 a 4 . (3.10) 

has no positive multiple roots, then the steady state is locally asymptotically stable for some 
cr 6 [0, cr 0 ) and it is unstable for cr > cr, with cr 0 < cr. At cr = cr 0 the Hopf bifurcation 
occurs. 


(ii) If 


{ 2 / 

a > ab min j 1, — 


(3.11) 


or 


a < Q'Zimin|l, and (a 2 + 2a 2 b 2 ^ < 2l(2aby + a(a 2 b 2 -/? 2 )) , (3.12) 

then at most one stability switch of the steady state is possible. Moreover, if the steady state 
is locally asymptotically stable for cr - 0, then it is stable for some cr e [0, cr 0 ) and it is 
unstable for cr > cr 0 , For cr = cr 0 the Hopf bifurcation is observed. 


Proof : Let m M - max{mi, m 2 }. Then the characteristic function (3.6) has the following form 


W(A) = 


1 


A (a + A)" ,M + 


(a + A) mM 

+ |/l(ra mi (a + A) mM ~ m + aba" 12 (a + A)'” 2 ) + ya m fa + A) 


rn M ~m\ \ - Act 


The zeros of W ( A) are the same as the zeros of 

D(A) = A 2 (a + A)' nM + ( A(ra mi (a + A)'"^'” 1 + aba m2 (a + A) mM “'" 2 ) + ya’ n< (a + A)'" M ~ m ' j e~ A<T . (3.13) 


For A = im define an auxiliary function 

F(oj) = oj 4 {a 2 + co 2 ) mM - || iaj(ra mi (a + ioj) mM ~ m ' + aba m (a + itof") + ya m '(a + ioj) mM ~ m 'f . 

(3.14) 

In the following we show, with one exception, that there exists a unique single positive zero = sfuf 
of the function F such that F\co 0 ) > 0. This, together with Theorem 1 from (45| allows us to deduce 
that the zeros of D crosses imaginary axes from left to right with a positive velocity. This yields the 
existence of stability switches and the occurrence of the Hopf bifurcation in the appropriate cases. 
For m 1 = m 2 = m, and for u = co 2 the auxiliary function F takes the form 


F(u) = u 2 (a 2 + uj - a lm (uj3 2 + y 2 ). 


(3.15) 


We are interested in the existence of the real positive roots of ( |3. 15| ). Because the number of sign 
changes between consecutive non-zero coefficients is equal to 1 the Descartes’ rule of signs indicates 
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that the polynomial ( |3.15[ ) has one simple real positive root denoted by uq. Clearly, the fact that 
coefficient of u 2+m is positive and F( 0) < 0, implies F'{uq) > 0. 

Moreover, 

F\u) = 2 u[a 2 + uj + mu 2 (a 2 + w ) m_1 - a lm /3 2 

and thus, F"(u) > 0. Hence, F'(«o) > 0. This completes the proof of part (a). 

For m\ = 1, m 2 = 2 and u = or the function F given by {] 3.14 [) has the form 


F{u) = u 4 + 2 a 2 u 3 + a 2 (a 2 - r 2 ^u 2 - a 2 (a 2 ft 2 + 2 a aby + y 2 ) u- 


a 4 y 2 


(3.16) 


We show that the coefficient of u in (|3. 1 6[) is negative and hence the Descartes’ rule of signs implies 


that a single stability switch for the trivial steady state of ( |3.1[ ) occurs. This is equivalent to the 
inequality 

G(a ) = a 2 /3 2 + 2 a aby + y 2 > 0, 

for a > 0. For b > n (that is the condition required for the existence of the trivial steady state of ( |3.1[ )) 
the discriminant of G is the following 

4y 2 ( a 2 b 2 - /? 2 ) = 4y 2 [a 2 b 2 - (r + ab ) 2 ) = -4ry 2 (r + lab) < 0. 

Thus, G(u) > 0, so F has exactly one simple positive root u 0 . Moreover, F( 0) < 0 hence F'(u {) ) > 0 
and the proof of part (b) is completed. 

Consider the case: m\ =2 and m 2 = 1. Then for u = to 2 Eq. ( |3.14[ ) takes the form ( |3.10| ) First, 
note that for the auxiliary function given by ( |3.10| ) inequality F{ 0) < 0, holds. Thus, F has at least 
one real positive zero. Hence, if F has only simple positive zero, by Theorem 1 from [[45]] the steady 


state of system ( |3. 1 [ ) is unstable for sufficiently large cr, which completes the proof of part (c.i). 

Now, we prove statement (c.ii). Assume that ( |3.11[ ) holds. Then, if a > ab, then the coefficient 
of u 2 is non-negative and independently on the sign of the coefficient of u the Descartes’ rule of signs 
indicates that there is exactly one simple positive real zero of F(u). Thus, at most one stability switch 
of the steady state can occur. On the other hand, if condition a > 2ayb/j3 2 holds, then the coefficient 
of u in F(u) is non-positive. Hence, independently of the sign of the coefficient of u 2 the single change 
of the sign is observed. Thus, F(u) has exactly one simple real positive zero. 


Now assume that (3.12) holds. To shorten notation denote 


a | = la by - a/3 2 


a 2 = a 2 b 2 


The first inequality of ( |3.12[ ) is equivalent to a\ >0 and a 2 > 0. Denote F^C) 
£ = u/a. Clearly, F(u) = 0 is equivalent to F\(() = 0. The function F i reads 

F\(Q = £ 4 + 2 < 3^ 3 - a 2 £ 2 + gai — y 2 . 


(3.17) 
F(a£)/a 4 , with 


We show that under the assumption ( ]3. 12[ ) the function F { is a strictly monotonic function of £, which 
implies that F is a strictly monotonic function of u. This, together with the fact /*’(()) < 0 implies that 
F has exactly one simple positive zero so the assertion of the point (c.ii) is true. 

In order to show monotonicity of F \, we prove that its first derivative is positive. We have 

F[(0 = 4? + 6a?-2a 2 £ + ai. 

The assumption a\ > 0 implies F\ (0) > 0. Now, we derive a condition that guarantees the positivity of 
F\ for all ^ >0. To this end, we calculate the minimum of F\(f) on the interval [0, +oo). Calculating 
the second derivative of F\ we find that the minimum of F[ is reached at the point 


C = + ^ a2 + foi)- 
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Using the fact that F"(£) = 0 we have 


F[(0 = ai - 2^ 3 - a 2 £. 

Thus, Fi is strictly increasing for all £ > 0 if and only if 

ill 2 + a 2 )f<a 1 . 


(3.18) 


A few algebraic manipulations indicate that ( |3 .1 8[ ) is equivalent to 

\ 3/2 


/ 2 2 V ( 2 ) 

I or + -a 2 I - a(a + ar 2 J<ari. 


Using the definition of a 2 the above inequality reads 

—^—(a 1 + 2a 2 b 2 ) 1 <a\+aa 2 b 2 . 

3V3 V ’ 

Due to the assumption a\ > 0 both sides of the above inequality are positive, so squaring that inequal¬ 
ity we obtain 

(a 2 + 2 crZ? 2 ) < 2l(ai + a a 2 b 2 ^ . 

Now using the definition of aq we get 

(a 2 + la 2 b 2 J < 2l{2aby + a(a 2 b 2 -/3 2 )) , 

which is fulfilled due to the second inequality of (|3.12[). This completes the proof of part (c). ■ 


Note, that there exist a set of parameters of system (3.1) with the shifted Erlang probability den¬ 


sities given by (3.4)—(3.5 > with m,\ = 2, m 2 = 1 such that the steady state is locally asymptotically 


stable for cr = 0 and conditions (3.12) hold. As a proper example we formulate the remark below. 


Remark 3.6 If 


2y 

ab < yS, and 1 < —, 


(3.19) 


then for all a > kfi+yr-cab the steady state of system (|3.1[) with the shifted Erlang probability densities 


given by ( |3.4[ )-( [33| ) with m\ = 2, m 2 = 1 loses its stability at cr = cr and the Hopf bifurcation occurs 
at this point. 

Proof : First, note that under the assumptions of the remark if a > ab min {1, ^], then condition 
( |3.1 1| ) holds and Theorem 


3.5 


implies the assertion of the remark. If a > + -j - ab and conditions 

(|3 .1 2[) hold then Theorem 33 also implies the assertion of the remark. We show that for some set of 


parameters these conditions can be fulfilled simultaneously. 


Consider a < aFmin jl, ^rj. The second inequality of Q3. 19[ ) implies that min jl, = 1. Hence, 


the stability condition of steady state for cr = 0 and the first condition Q3. 12[) reads 


1 2y 

-B -1 - ab < a < ab. 

2 h /3 


(3.20) 


First, we show for some parameters inequalities Q3.20[ ) give a non empty set of a. This is true if the 
following inequalities 


1 2y 

—B H-< 2 ab < 2J3 

2 h /3 


(3.21) 
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hold, where the second inequality in ( |3.21[ ) is just the first one in ( |3.19| ). Inequalities ( [3.21[ ) do not 
contradict each other if and only if 

1 2 y 

-/3+^-<2/3. (3.22) 


Thus, if Jr < 3/2, then there is a non-empty set of ab such that inequalities ( ] 3.21 [ ) hold and for such 


ab (3.20) determines a non-empty set of a. 


Now, we show that if (3.19) and the first inequality of (3.12) hold then the second inequality 


of (3.12) also holds. 


Note that the assumption ab < (3 means that the right hand side of the second inequality of ( |3.12| ) is 
a decreasing function of a (of course we need to use here the first inequality of ( |3.12| ) that guarantees 
the positivity of the expression under the second power), while the left hand-side is an increasing 


function of a. Since our assumptions imply that a < ab, it is enough to check if inequality ( |3.12[ ) 
holds for a = ab. 

Rearranging terms we obtain 

27 a 6 b 6 < 2l{o?b 3 + ab(2y -/3 2 jf. 

The above inequality obviously holds as 2y > ft 2 , which completes the proof. ■ 

For the general case of system ( |3.1[ ) with the shifted Erlang probability densities given by ( ]3.9[ )~ 


( |3.5[ ) the characteristic function D(A) is given by ( |3.13[ ) and the auxiliary function F is given by ( |3.14[ ). 
It is easy to see that the highest power of co is 4 + 2 m M and the coefficient of it is 1. Moreover, we 
have F{ 0) < 0. Thus, there exists oj () > 0 such that F(co 0 ) = 0 and the roots cross the imaginary axis 


from the left to the right half-plane. Thus, using Theorem 1 from [45], we can deduce the following 


Remark 3.7 If all roots of F((o) given by ( |3.14[ ) are simple, and the trivial steady state of system ( |3.1[ ) 
with the shifted Erlang probability densities given by — ( |3.5| ) is locally asymptotically stable for 
cr - 0, then it loses its stability due to the Hopf bifurcation for some cr 0 > 0 and it is unstable for 
cr > era, with some cr > cr 0 . 

It seems that the case of multiple roots of F(oj) is non-generic, however we will not prove it here. 


3.2.2 Piecewise linear probability densities 


For the function defined by (3.3) we have 


poo <x 

1 ^-**=5* 


(e' )£ + e -' l£ 


2 £-' lcr 

- 2 ) = -^r( c °sh(/le) - l) 


and thus the characteristic function for the trivial steady state of system ( |3.1[ ) reads 

X oo s~*oo 

(r/i(r) + abf 2 {T))e~ AT dr + y J f(j) e"' lT dr, 

which is equivalent to 

2 ~-<icr 2 p~' lcr 

VT(T) = A 2 + (Ar + y) ,, (cosh(/lg) - l) + Aab —— r (cosh(/le) - l). 


(3.23) 


Finding zeros of function ( |3.23[ ) is not a trivial task. If the characteristic function has the form W(A) = 
P(A) + Q(A) c~ Ar , where P and Q are polynomials, one could use the Mikhailov Criterion [46], to 


estimate the number of zeros of the characteristic function lying in the right half of complex plane. 
However, in the considered case we do not have strict polynomials. Therefore, we use a generalised 
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Mikhailov Criterion (for details see A). Clearly, cosh(z) = o > for z e C is an analytic function, 
hence cosh V ) ^ 1 is also analytic. Consequently, the function W defined by (|3.23[) is also analytical. 


Moreover, obviously condition ( |A.1[ ) holds. Thus, we can apply the generalized Mikhailov Criterion. 
Considering cr > s > 0, we have 


W(X) = A 2 + {Ap + y)g{Ae)z Acr , 


(3.24) 


where 


Thus, 


2 (cosh x - 1 ) , ^ , 2 - 2 cos* 

g(x) =- 5 -, gi(x) = g(ix) =- - -. (3.25) 


Re(W(zco)) = -to 2 + gi(to£)(y cos(tocr) + poo sin(cocr)). 


(3.26) 


Im(VT(za>)) = gi(tos)(j3to cos(tocr) - y sin(cucr)). (3.27) 

Here we formulate sufficient condition for stability of the trivial steady state of system ( |3.1[ ) for 
the piecewise linear distributions defined by ( |3.3[ ) with cr > s > 0. This condition is clearly not 
necessary, what we will show numerically in the next section. 


Proposition 3.8 Let b > p, p and y be defined by ( |3.7| ). Then for 

(i) 


cr < 


n 


(3.28) 


P + y/p 2 + 4y + |7T 

the trivial steady state of system ( |3. 1 [ ) with the piecewise linear distributions defined by ( ] 3.3 [ ) 
and cr > s > 0 is locally asymptotically stable. 


(ii) 


/? In 

- <cr < - 

y fi+ y/fi 2 + 4y 


(3.29) 


the trivial steady state of system ( |3.1| ) with the piecewise linear distributions defined by ( |3.3[ ) 
with cr > s > 0 is unstable. 

Proof : We use the generalized Mikhialov Criterion (for details see |A]) to show that the change of 
the argument of W(ico) as to varies from 0 to oo is equal to n. In fact, we show that the behaviour of 
the hodograph is more restrictive, namely that for some do > 0 we have Re(VT(/o»)) < 0 for oj > do, 
while Im (W{ico)) > 0 for to e (0, do\. To this end first, note that if gftoe) = 0 and oo > 0, then 
Re(Wf7m)) < 0. Thus, without loss of generality we assume that gfcoe) > 0. 


First, let us estimate the real part of W(ito) given by p.26[ ). Since b > p and gi(ecu) < 1 we have 

Re(VT(zcu)) < -to 2 + y + top =: gR{to). 

Now, let us consider the imaginary part of W(ito). In fact, since we are interested in the sign of it and 
gi(ecu) > 0 it is enough to consider the sign of 


where 


top cos (tocr) - y sin(cucr) = gfitocr), 


gi(x) = — x cos x — y sin x. 
cr 


Note, that ( |3.28| ) implies p > ycr, and thus gi(x) > 0 for x close to 0. On the other hand, gi(n/2) = 
-y < 0. It is also easy to check, that g/(x) has a unique zero in (0, n/2) for p > ycr. Below, we give 
some estimation for this zero. 
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Since tan jc < \ -rr x , for x £ (0, n/2) we have 


B n x 

—x > --> tanx. 

ycr 2%-x 


(3.30) 


A simple algebraic yields that the first inequality of (3.30) is fulfilled for all 


n 


0 <x< -|1 


CT • 


< 


n 


Thus, as long as cocr < | (l - cr • we have Im(lT(uu)) > 0. 

However, g R is a quadratic polynomial with g#(0) > 0. Moreover, the coefficient of or is negative 
and the positive root of g R is equal to \{j3 + \jfi 2 + 4y). An easy algebraic manipulation shows that 
condition (|3.28[) is equivalent to 




+ 


V^ + 4y)<^(l- 


r 

a ' P 


and thus g R {co) < 0 for all crco > ?(l - cr • |j, which completes the proof of part (i). 

The idea of the proof of part (ii) is similar to the proof of part (i). However, this time we want 
to show that for era* 6 (0, n) the imaginary part Im (W(iuSj) < 0, while for crco > n the real part 
Re(VTO'm)) < 0. This indicates that the hodograph goes below the (0,0) point on the complex plane 
when co tends to +oo and hence the change of the argument of W(ico) differs from n. Hence, from the 


generalized Mikhialov Criterion it is clear that the trivial steady state of system ( |3.1[ ) is unstable. 

Note, that for cocr £ {n/2, n) we have PcCWdoo)) < 0 due to the fact that on this interval cosine 
is negative and sine is positive. Consider cocr £ (0, ;r/2). Thus, cosine is positive and inequality 
fico cost cocr) - y sin(mcr) < 0 is equivalent to 


—cocr < tan(cucr). 
cry 


(3.31) 


The right hand side of (331) is a strictly increasing convex function on (0, n/2) and have the first 


derivative at 0 equal to 1. Thus, the first inequality of ( |3.29[ ) implies p.31[ ). 

Now, it is enough to show that g R (n/cr) < 0. A straightforth calculation yields 


l n \- n 

- - +fi— + y < o. 
\cr / cr 


Solving quadratic equation we get 


cr < 


2 n 


(3+ y/f3 2 + 4y ’ 


which is the second inequality of (3.29). 


Remark 3.9 Note, that ( |3.28[ ) implies cr < ply and consequently, conditions ( |3.29[ ) and ( |3.28[ ) are 
mutually exclusive. 


Theorem 3.10 If the trivial steady state of system (3.1) with the piecewise linear distributions defined 


by ( ]3.3[ ) and a > s > 0 is locally asymptotically stable for cr - s and the function 

F(co) = co 4 - s 2 (arp 2 + y 2 )g 2 fsco), 


where g\ is defined by ( ] 3.25 [ ) has no multiple positive zeros, then the steady state is locally asymptot¬ 
ically stable for some cr £ [<s, <x 0 ) and it is unstable for cr > cr, with cr 0 < cr. At cr = cr 0 the Hopf 
bifurcation occurs. 
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Proof : For the considered case the characteristic function is given by ( |3.24[ ). It can be easily seen 
that if W{ioj) = 0, then 


- ito 1| 2 = e 2 ||hn/? + y||~ g(ica £) 


and thus 


a> 


= s 2 (u 2 /3 2 + y^ffsoj ), 


where is defined by ( ]3.25| ). The function F has the following properties 

F( 0) = -s 2 y 2 < 0, and lim F(oj) = +oo. 

CO —> + 00 

Hence, we conclude that there exists a > 0 > 0 such that F(co 0 ) = 0 and F'(co 0 ) > 0. It is assumed that 
all zeros of F(a>) are simple, thus F'(a> o) > 0 and by Theorem 1 from [451, the thesis of theorem is 
proved. ■ 


Remark 3.11 Note, that the case where the function F has multiple positive zeros is not generic. 

Proof : Consider the case of the multiple positive zero(s) of F. First, note the following facts: 

(a) F is a analytic function so it has isolated zeros; 

(b) For (o > maxfe a/# 2 + y 2 ,1} inequality F(a>) > 0 holds. 

(c) If gi(eo») = 0, then F(to) = ( Ikn/s) 4 > 0 for some k e N, k > 1. 

Facts (a) and (b) imply, that F has a finite number of zeros. Suppose that for the arbitrary values of fio 
and y 0 the function F has the multiple zero at (o m j, for some j = 1,2, ..Jm, where j M > 1 is a natural 
number. Now, consider the interval / = [0, maxje y//3 2 + y 2 ,1} + 1]. Clearly, dF/dy < 0 holds for all 
oj ■+■ Ikn/s. Moreover, F' is also an analytic function so it has a finite number of zeros in the interval 
/. Thus, we choose 0 < e y < 1 so small that for all y 6 (y 0 - e r , y + e r ) \ {y 0 } the function F has no 
multiple root in the interval /. This completes the proof. ■ 


4 Stability results for parameters estimated by Hahnfeldt et al. 


In the previous section, we analytically investigated the stability of the trivial steady state of the family 
of angiogenesis models with distributed delays p.l[ ) and distributions characterized by the probability 
densities f given by ( ]3.3[ ) and ( ]3.4[ )-( |33| ). System ( |3. 1[ ) is a rescaled version of ( ]2. 1 [ ), where the trivial 
steady state of ( |3.1[ ) corresponds to the positive steady state of ( |2.1[ ). 

In this section we illustrate our outcome for particular model parameters considered earlier in the 
literature, that is the parameters estimated by Hahnfeldt et al. in the case of the model without the 
treatment (see [47;]). Namely, we consider the following set of parameters: 


p = 0, a H = 8.73 x 1(T 3 , b = 5.85, r = 0.192, 


(4.1) 


and we take the function h{9) = -In 9. We present stability results for the distributed models for 


two different values of parameter a. To keep the description short, for model p.l[ ) with a = 1 we 
would refer to as the Hahnfeldt et al. model, while to model (13. lb with a = 0 as the d’Onofrio- 


Gandolfi model. However, in this paper we also consider positive values of p that can be interpreted 
as a constant anti-angiogenic treatment. 
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Figure 1: Dependence of the critical average delay value in the case of the non-shifted Erlang distribu¬ 
tions given by ( |3.4[ )-( |33| ) with ni\ = m 2 = m on the constant treatment coefficient /a for system ( |3. 1 [ ). 
The critical average delay is defined by r cr ,o = m/a CIt o, where a crj0 is the critical value of parame¬ 
ter a for which stability change occurs (see Proposition |3.3[ ). In the left-hand side panel results for 
the Hahnletdt et al. model (a = 1) are shown while in the right-hand side panel results for the 
d’Onofrion-Gandolfi model (a = 0). All other parameters, except /j. are defined by (|4. 1 [). 


4.1 Erlang probability densities 


First, focus on model p.l| ) with the non-shifted Erlang distributions given by ( |3.4| )-( |33] ) with cr = 0, 
i = 1,2. The expression m/a describes the average delay. On the other hand, in the more general case 
of m\ 4 m 2 a similar interpretation leads to the conclusion that the average delays are mi/a and m 2 /a. 
One could also consider the arbitrary values a = a\ for f\ and a = a 2 for / 2 , but then the analytical 
expressions become very complicated and we decided not to study these cases here. 

In Fig. |T] the dependence of the critical average delay value r cr 0 with m, = in (i = 1,2) on the 
constant treatment coefficient /a is presented. In the case of the non-shifted Erlang distributions the 
average delay is defined by r cl! o = m/a cr ,o an d the left-hand side panel and the right-hand side panel 
present results for the Hahnfeldt et al. model (a = 1) and the d’Onofrion-Gandolfi model (a = 0), 
respectively. The curves were calculated from the stability conditions given in Proposition |3.3[ The 
stability regions of the trivial steady state of ( |3.1| ) with the non-shifted Erlang distributions are below 
curves, while above them there are the instability regions. It is worth to emphasise that there is 
a qualitative difference between the cases: m = 1 and m = 2,3 for the Hahnfeldt et al. model, 
while there is no such a difference for the d’Onofrio-Gandolfi model. Clearly, for the Hahnfeldt et 
al. model and m = 1 the critical average delay r cr o tends to +oo as /a —* b = 5.85 while for m = 2 
and m = 3 we have r cr , o\ b = 4//? « 0.662 and r^o) fc = 8 /(3/3) « 0.441, respectively. Nevertheless, 
for the d’Onofrio-Gandolfi model the critical average value of delays are almost the same for the case 
m =1,2 and m = 3 for // 6 [0, b) and the behaviour is qualitatively comparable with the result for the 
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Figure 2: Dependence of the critical value of parameter a for the trivial steady state of system ( |3.1| ) 
with the non-shifted Erlang probability densities given by ( |3.4[ )-( |33| ) (i.e. cr = 0). The regions above 
the plotted curves corresponds to the stability regions for the particular choices of parameters m,, 
i = 1,2. For m\ = 2, m 2 = 1 and a = 1 the trivial steady state is always stable. All other parameters, 
except /u, are defined by ( |4.1[ ). 


Hahnfeldt et al. model for m = 1 .Moreover, in the considered cases the critical average delay values 
in the case of the non-shifted Erlang distributions are the increasing functions of /u. This means that 
the increase of treatment increases the stability region. 

In Fig. [2] we plot the critical value of parameter a, that is a CTt o, for the trivial steady state of sys¬ 
tem p.l[ ) with the non-shifted Erlang probability densities given by ( ]3.4| )-( [33| ) (i.e. cr = 0) calculated 
based on Proposition |3.3| Here, the regions above the plotted curves correspond to the stability regions 
for the particular choices of parameters m,, i = 1 , 2 , while those below correspond to the instability re¬ 
gions. It should be clarified that for the case a - \,m\ =2 and m 2 = 1 the steady state is stable for all 
values of a . We see that all plotted functions are decreasing as functions of the strength of constant 
treatment [j. However, already for the non-shifted Erlang probability densities we observe a large 
differences between both models regarding the model dynamics. For the Hahnfeldt et al. model we 
have the largest stability region for m { = 2 and m 2 = 1, while in the case of the d’Onofrion-Gandolfi 
model for m\ = 1 and m 2 = 2. Moreover, for the case ni\ = 2 and in 2 = 1 the steady state is stable 
independently of the parameter a for the Hahnfeldt et al. model, while in the same case of m, for the 
d’Onofrio-Gandolfi model the steady state might change its stability with the change of a. Addition¬ 
ally, for a = 0 the sizes of the stability regions for m, = 2 and m\ = 2, m 2 = 1 are the same and it is 
not the case a = 1 . 

The stability of the steady state for the non-shifted Erlang distributions can be determined by the 
Routh-Hurwitz criterion, although it may require tedious calculations, in particular for large values of 
parameters m\ and/or m 2 . On the other hand, the analytical results we obtained for the shifted Erlang 
distributions are rather limited. The results presented in Theorem [33] are only the existence results and 
give no information of the magnitude of the critical value of cr for which stability is lost. Although, it 
is possible to derive an expression for the critical value of cr, it would involve arccos of some algebraic 
function of a>o, which value cannot be, in general, determined analytically and the final result would 
not be informative. In consequence, we decided to calculated numerically the critical values of the 
average delay r c;vr for the considered model with the shifted Erlang distributions for parameters given 
by ( ]4.1[ ). Clearly, if one wish to calculate the critical average delay for model with the shifted Erlang 
distributions it appears that it depends on sigma directly and is given by t C[(T = m/a + cr cr , where 
cr cr is a critical value of cr defined in Theorem 3.5 In Fig. [3] we present the stability results for the 
trivial steady state of system ( |3.1[ ) with the shifted Erlang probability densities given by ( |3 .4[ )— ( f!T75| ). 
In particular, in Fig. [3] (left column) we plot the dependence of r c „ (r on the treatment strength /a for 
the particular choices of parameters m ; and the particular choices of parameter a for the model with 
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Figure 3: Dependence of the critical average delay t cjjt for the trivial steady state of system ( |3. 1 [ ) 
with shifted Erlang probability densities given by ( |3.4| )-( |33] ) for the particular choice of parameter a 
and comparison with the critical value of the average delay r cr> 0 in the case of cr = 0. For the case of 
the Hahnfeldt et al. model (i.e. a = 1), presented in the left column, the regions above the plotted 
curves correspond to the instability regions for the particular choices of parameters m\ = m 2 and of 
parameter a. The stability region is the region below curves. For the d’Onofrio-Gandolfi model (i.e. 
a = 0) presented in the right column, the difference between critical value of the average delay r cr 0 
in the case of cr = 0 and the critical average delays r cr/r for chosen values of a are so small that the 
plot would be non-informative. Thus, in the right column we presented the difference cr+ = T cl - () -T cr (r . 
All other parameters, except /r, are defined by ( |4.1| ). 
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a = 1. Similarly as in Fig. |T| the regions below the plotted curves correspond to the stability regions. 
Additionally, we plot T cr fi(p) curves to indicate thresholds for which the destabilizations occur in the 
case cr = 0. If the curve for a considered a is above the T rr (l curve, then the steady state is unstable for 
cr = 0 and remains unstable for all cr > 0. In the both cases m, = 1 and m, = 2 the increase of the value 
of the parameter a (which is equivalent to decrease of the critical average delay), defining the shape 
of the probability densities, implies the decrease of the stability regions of the steady state. Since 
the corresponding curves for cr = 0 are very close to each other instead of plotting them directly we 
decided to plot the differences cr + = r c ,, 0 - r clvr between the critical values of the average delays r c „ 0 
in the case cr = 0 and the critical average delays r cr cr for the chosen values of a. Clearly, whenever 
plotted curve reaches the jj axis the trivial steady state becomes unstable for all cr > 0. From Fig. [3] 
(right column) we deduce that the stability areas again decrease with the increase of the parameter a. 


4.2 Piecewise linear distributions 


Let us focus on the stability and instability regions of the trivial steady state of system (J3TT|) with 
piecewise linear distributions defined by ( |3.3[ ). Since in such a case we have smaller number of 
parameters defining the shape of the probability density we are able to plot the stability regions in 
(cr, s ) plane, see Fig. |4j Note that, for s > cr system ( |3. 1[ ) becomes a neutral system, which is out 
of our consideration hence this region is greyed. Clearly, the estimations obtained in Proposition |3.8| 
are rough. However, solving numerically a system ReVF(zaf) = 0, ImVFf/oj) = 0, where real and 
imaginary part of characteristic function are given by ( |3.26[ ) and ( |3.27[ ), respectively, we are able to 
calculate the stability region and the curve for which the stability change occurs. This line was plotted 
in Fig. [4] as a thick solid blue line. 

Additionally, we plot (in both panels) the vertical dashed lines that denote the conditions guar¬ 
anteeing stability (see Proposition |3.8( i)) or instability (see Proposition |3.8[ ii)) of the trivial steady 
state. For the Hahnfeldt et al. model (left panel in Fig. [4]) the condition ( |3.28[ ) from Proposition |3.8| 
is denoted by a vertical dashed line. In that case, the condition (|3.29[) cannot be fulfilled since the 


expression on the right-hand side of (3.29) is always smaller than (3/y. For the d’Onofrio-Gandolfi 


model the dashed vertical line on the left to the solid vertical line indicates the condition (3.28). Two 


dashed vertical lines on the right to the solid one indicate the region in which the condition ( |3.29[ ) is 
fulfilled. Those numerical results show that the conditions from Proposition |3.8| are only sufficient 
but not necessary. From the figure we also deduce that the stability region for the trivial steady state 
for the d’Onofrio-Gandolfi model is smaller than the one for the distributed Hahnfeldt et al. model. 
Moreover, computed by us the critical values of cr for s —> 0 and both a = 1 and a - 0 agrees with 


those calculated by Piotrowska and Forys in [33J for the models with discrete equal delays. That 
agrees with the intuition, since for e —> 0 models given by (|3.1[) with piecewise linear distributions 


defined by (3.3) reduce to the models with double discrete delay. 


For /j > 0 and the Hahnfeldt et al. model we observe a move of the stability switch curve to the 
right with hardly any change in the shape indicating that the stability region increases with increase 
of parameter /u, see left panel in Fig. [5] This shift is very small. Moreover, for /j = b we obtain 
limiting values: 0.26 for s —> 0 and 0.33 for s —> cr. On the other hand, for the d’Onofrio-Gandolfi 
model the change is more visible. The lines start to lean to the right. For small //, the change is quite 
small, compare with right panel in Fig. [4j For example, for /r = 3 we obtain cr = 0.509 for s —> 0 
and cr = 0.59 for s = cr. As /r approaches to b the changes become more rapid. For /r = 5.7 we have 
cr = 5.326 for s —» 0 and cr = 5.632 for e - cr, while for n - b we have cr = 8.181 and cr = 10.065 
for e —» 0 and e-ct, respectively. 
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Figure 4: Stability and instability of the trivial steady state of system ( |3.1| ) with piecewise linear 
distributions defined by ( ]3.3[ ) with cr > 0, s > 0, i = 1,2 for the parameters given by ( [4. 1 [ ) with 
a = 1 (the left hand-side panel) and a - 0 (the right-hand side panel). The greyed part denotes the 
region where s > cr for which system ( |3. 1 [ ) becomes a neutral system. The solid lines indicate the 
stability switch. In the left-hand side panel, the vertical dashed line denotes the condition ( |3.28| ) from 
Proposition |3.8| In the right-hand side panel, the vertical dash line in the stability region denotes the 
scope of the condition p.28[ ) form Proposition |3.8[ while two vertical dash lines in the instability 
region denote the scope of condition ( |3.29[ ). 




Figure 5: Stability and instability of the trivial steady state of system ( |3 .1 1 ) with piecewise linear 
distributions defined by ( |3.3| ) with cr > 0, s > 0, i = 1,2 for the parameters given by ( |4.1| ) with 
a = 1 (the left hand-side panel) and a = 0 (the right-hand side panel) and different /i. The greyed 
part denotes the region where s > cr for which system ( ]3. 1[ ) becomes a neutral system. The solid lines 
indicate the stability switch for (from left to right) /a = 0,3, b for a = 1 and for // = 0,3,5.7, b for 
a = 0. 


5 Discussion 


Although discrete delays often appear in models describing various biological phenomena, e.g. [48 


58] and references therein, it is obvious that in natural processes delay, if it exists, is usually somehow 
distributed around some average value due to differences between individuals and/or environmental 
noise. Thus, we believe that distributed delays are more realistic. However, in some cases, when 
the distribution has small variance and compact support discrete delays may be treated as a good 
approximation. Clearly, models with discrete and distributed delays might have different dynamics in 
some ranges of parameters. In the presented paper we have compared the behaviour of two types of 
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systems in the context of the stability of the steady state for a particular family of models describing 
the tumour angiogenesis process. 


Table 1: The critical average delay value r cv .o in the case of non-shifted Erlang distributions given 
by ( |3.4[ )-( |33T ) with my = m 2 = m and a\ = a 2 = a for arbitrary chosen parameters n for system ( |3.1[ ). 
The critical average delay is defined by r cl2l = m/a c r ,o, where a cr ,o is the critical value of parameter a 


for which the stability change occurs (see Theorem [33]). 



Hahnfeldt et al. model (or = 1) 

d’Onofrion-Gandolfi model (a = 0) 

b 

0 

2 

4 

5.85 

0 

2 

4 

5.85 

Illy = m 2 = 1 

8.069 

12.261 

25.515 

OO 

0.256 

0.390 

0.811 

OO 

niy = m 2 = 2 

0.611 

0.628 

0.645 

0.662 

0.253 

0.382 

0.780 

20.833 

niy = m 2 = 3 

0.421 

0.428 

0.435 

0.441 

0.252 

0.380 

0.770 

13.889 


Presented analysis of distributed models includes: the uniqueness, the positivity and the global 
existence of the solutions, the existence of the steady state and the possibility of the existence of the 
stability switches. We have analytically derived conditions, involving the parameters defining the 
probability densities, guaranteeing the stability or instability of the steady state. We have also shown 
that in some cases the single stability switch is observed and the Hopf bifurcation occurs. 

For the particular set of parameters, estimated by Hahnfeldt et al., we have investigated the sta¬ 
bility regions for steady state. We compared the results for both the Hahnfeldt et al. (a = 1) and 
the d’Onofrio-Gandolfi (a = 0) models in the case of different probability densities. For both models 
we considered the Erlang shifted and non-shifted distributions as well as the piecewise linear dis¬ 
tributions. We want to emphasise here that, in the general case, it is hard to say for which model 
Hahnfeldt et al. or d’Onofrio-Gandolfi the stability region is larger since it strongly depends on the 
considered probability densities and theirs shapes. However, we observe certain similarities. 

First, for /a = 0, i.e. the family of models without treatment, we see that for the Hahnfeldt et al. 
model with non-shifted Erlang probability density the larger m, = m are the smaller the stability region 
is, see Fig. [I] and Table [T] The same holds for the d’Onofrio- Gandolfi model with the non-shifted 
Erlang distributions, see Tableland moreover, for all considered m, = m, i = 1,2 (and /a = 0) for the 
Hahnfeldt et al. model the stability region is larger than for d’Onofrion-Gandolfi one. Similarly, for 
the models with the piecewise linear distributions the stability region for /i = 0 for the Hahnfeldt et 
al. model is larger than the one for the d’Onofrio-Gandolfi model, compare Fig. [4] 

If we consider a positive parameter p smaller than b (to ensure the non-negativity of the steady 
state of model ( |2.1[ ), which corresponds to the trivial steady state of ( |3.![ )), we observe a similar 
tendency for the Erlang distributions in dependence on the parameters niy = m 2 (compare Table jTj 
for the non-shifted distribution case). However, the dependence on /j depends strongly on the chosen 
model. For the Hahnfeldt et al. model and my = m 2 > 2, this dependence is almost linear and very 
weak, while for the d’Onofrio-Gandolfi model it is much stronger. In results, for my = m 2 > 2 and any 
sufficiently large value of n e [0, b ) the critical average delay for the Hahnfeldt et al. model becomes 
smaller than for the d’Onofrio-Gandolfi model. The case niy = m 2 = 1 is different. Dependence on 
H is similar for both models and the critical average delay for the Hahnfeldt et al. model stays larger 
than for the d’Onofrio-Gandolfi model for any given value of /a. Similarly, for different values of m ( , 
it seems that dependence on /u is stronger for the d’Onofrio-Gandolfi model, see Fig. [2] Nevertheless, 
for the d’Onofrio-Gandolfi model with the non-shifted Erlang distributions there is no difference 
regarding the stability results between cases niy = m 2 = 2 and my = 2 , m 2 = 1, which is not the 
case for the Hahnfeldt et al. model, where for niy = 2, m 2 = 1 we have stability independently on 
the value of parameter a. For the shifted Erlang distributions we have investigated the changes of the 
size of the stability region also in the context of the change of the value of parameter a. For both the 
Hahnfeldt et al. and d’Onofrion-Gandolfi models we see that for all considered m = m ; , i = 1,2, the 
increase of the parameter a decrease the stability region, however in the case of d’Onofrio-Gandolfi 
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model we have compared the differences of the T cr ,o - r cr iT . A similar increase of the stability region 
with a decrease of concentration of delay distribution was observed in [ 361 for one linear equation 
with distributed delay and Erlang probability density. 

Nevertheless, we see that since for all considered cases of the shifted and non-shifted Erlang dis¬ 
tributed models t civt and T cr {h respectively, are increasing (sometimes slightly) functions of variable /a 
implying that the increase of the constant treatment strength enlarges the stability regions for the pos¬ 
itive steady state of model ( |2.1[ ). Our analysis also shows that the increase of the positive parameter 
H enlarges the stability area for the steady state for both (Hahnfeldt el al. and d’Onofrion-Gandolfi) 
distributed models with the piecewise linear distributions, but this time for the d’Onofrion-Gandolfi 
model this increase is more pronounced than for the distributed Hahnfeldt et al. model, see Fig. [5] 
Moreover, for small values of parameter /a the stability region for the Hahnfeldt et al. model with the 
piecewise linear distributions is larger than the one for the d’Onofrio-Gandolfi model, while for the 
larger values of fa the situation is opposite. Performed simulations show that the change occurs before 
/a * 1.42. 

The variances for the shifted and non-shifted Erlang distributions are given by and they give 
a measure of the degree of the concentration of the delay around the mean. Actually, a better measure 
of the spread of the distribution around the mean for our purposes is the coefficient of variation, i.e. 
the ratio of the standard deviation to the mean, that is y/mi/(acr + m,). Clearly, for the non-shifted 
distributions we have 1 / yjrn,, which implies that larger parameter m, yields smaller the considered 
ratio. Hence, the increase of m, decreases the percentage dispersion of the average delay. On the 
other hand, for the shifted Erlang distributions, the coefficient of variation is a decreasing function of 
cr. This dependence is obvious since the increase of cr with fixed m, and a means that the average 
delay is increased while standard deviation remains constant. The coefficient of variation is also 
a decreasing function of parameter a. On the other hand, if we fix a and cr, and study the influence of 
ntj, then the coefficient of variation increases if m, < acr and decreases otherwise. 

For piecewise linear distributions (cr > s) we have the average value of f (defined by ( |3.3[ )) equal 
to cr and the standard deviation given by e/ V6. Hence, the coefficient of variation equals to e/(cr V6). 
Thus, it is an increasing function of s (for fixed cr) and a decreasing function of cr (for fixed e). Thus, 
we conclude that the percentage dispersion of the average delay is the increasing function of s and 
decreasing function of cr and it is always smaller than 1. Clearly, all that should be taken into account 
whenever the probability distributions describing the characteristics of delays are estimated. 

We believe that with the proposed type of models one can describe in a more realistic way the 
process of angiogenesis. However, the considered model should be validated with the experimental 
data. Our results show that a system behaviour in the case of distributed delays might be different from 
the behaviour of the system for discrete delays. Hence, to validate the model with the experimental 
data it is important to choose the right type of the delay(s) in the model preferably based on the 
experimental data or hypothesis postulated by the experimentalists. Our study also shows that a shape 
of probability density for the distributed delays has an essential influence on the model dynamics. 
On the other hand, it is not a trivial task to choose the right distribution if one has not sufficiently 
large set of the experimental data and does not assume any particular shape of the distribution. This 
problem need to be individually addressed according the available data and it should be done for the 
considered model in the future. In our opinion, models with distributed delays could be also studied 
in the context of efficiency of different types of tumour therapies. 

We would like to emphasise that oscillations in the experimental data for different human and an¬ 
imal cell lines in the context of neoplastic diseases, even without administration of any treatment, has 
been already widely observed for leukaemia ( [59} 601), B-cell lymphoma ( pTTJ ) and solid tumours 
( [62}63|). To investigate the origin of such phenomenon a number of mathematical and numerical 
models has been developed. For example Kuznetsov et al., [641, studying the interactions between 
an immune system and a tumour by system of ODEs, pointed out that different local and global bi¬ 
furcations (observed for the realistic parameter values they estimated) showed that there might be 
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a connection between the phenomena of immunostimulation of tumour growth, formation of “dor¬ 
mant” state and “sneaking through” of tumour. Similarly, Kirschner and Panetta in [ 651 discussed 
inter alia the biological implications of the Hopf bifurcation in the case of the tumour-immune system 
interactions. More recently, Pujo-Menjouet et al. [ |66] , using the model with discrete delay for which 
the Hopf bifurcation is observed, explained why and how short cell cycle durations gave rise to long 
period oscillations of order form 40 to 80 day for periodic chronic myelogenous leukaemia patients. 
At the same time in [67) the bifurcation phenomenon and its role for a white-blood-cell production 
model with discrete delay was studied. The experimentally observed oscillations in the tumour vol¬ 
ume [ |62|[63| can be explained by the analytic results from earlier and more theoretical paper [ [68) 
focussing on the growth of avascular multicellular tumour spheroids, where the delay in the process 
of proliferation was considered. Later on, this work was extended in J69) , while the model with the 
delay only in the apoptosis process was considered in [70|. Next, the model that took into account 
delays in both processes was studied in [71]. It is believed that the p53 protein plays an essential 
role in preventing certain types of cancers and is called “guardian of the genome” J72) , hence we 
could also mention here the models of Hesl and p53 gene expressions, both with negative feedback, 
and with discrete time delays proposed by Monk [51]. The model of the Hesl protein was examined 
later by Bernard et al. [73] and Bodnar and Bartlomiejczyk [501.In that cases the mathematical anal¬ 
ysis showed, that experimentally observed oscillation can be explained by the time lag present in the 
system due to the DNA transcription process as well as the diffusion time of mRNA and the Hesl 
proteins. 

Finally, we discuss shortly possible generalisation of the results obtained in this paper. The first 
equation of the considered model is quite general and can describe a process with saturation. On 
the other hand, the second equation comes from a quasi-stationary approximation of some reaction- 
diffusion equations under particular assumptions that are justified for the angiogenesis process. Thus, 
most probably it cannot be straightforwardly transformed to describe other biological or medical pro¬ 
cesses. However, the stability analysis performed in the paper is based on the local properties of the 
functions. Thus, we may assume a more general form of the second equation namely q(t)G(p t /q t , p) 
(where p t and q, according to a functional notation typical for DDEs denote terms with delay) assum¬ 
ing that G is an increasing function in the first variable and decreasing in the second one. Then, in 
our opinion, under some additional assumptions, one should obtain results similar to those presented 
in the paper. 


A Generalized Mikhailov Criterion 


Here we formulate a generalized version of the Mikhailov Criterion (see eg. [46]). The classical 
formulation of the Mikhailov Criterion is for the characteristic functions that are sums of polynomials 
multiplied by an exponential function. However, it can be generalized for wider class of functions. 
Below, we present a detailed formulation. 


Theorem A.l (Generalised Mikhailov Criterion) Let us assume that W : C —» C is an analytic 
function, has no zeros on imaginary axis and fulfils 

W(A) = A n + o(A"), W'(A) = nA"- 1 + o(A n ~ 1 ), (A.l) 

for some natural number n. Then the number of zeros ofW in the right-hand complex half-plane is 
equal to n/2 - A/n where 

A = A wm+oo) argW(icL>). 

Note, that A denotes the change of argument of the vector W(ioj) in the positive direction of complex 
plane as oj increases from 0 to +oo. 
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Proof: The proof is exactly the same as the proof of the Mikhailov Criterion in [46, Th. 1], It is based 
on integration of the characteristic function on the following contour: part of imaginary axis, portion 
of the imaginary axis from ip to -ip and a semi-circle of radius p with the middle in zero located in 
the right half of complex plane. The fact that W is an analytic function implies that it has the isolated 
zeros and all integration can be done in the same way as in [ ]46| . Moreover, the assumption ( |A.1[ ) 
implies that all limits calculated in [46] remain the same. ■ 
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